#!/usr/bin/python -tt


######### standard library imports ###########################################
import struct
import sys
import time 
from pprint import pprint
import logging
import logging.config
import os
import rlcompleter
import readline
from optparse import OptionParser

######### external library imports ###########################################
import numpy as np
import matplotlib.pyplot as plt
from ROOT import *

######### project library imports ############################################
import new_helpers as helpers
from Turnable import *
import readcorsika
import absorb_photons
import reflector   
import fact_spline_interpolations as fact_si
import turnIntoZenit

scriptpath = os.path.dirname(os.path.realpath(__file__))
logging.config.fileConfig(scriptpath+'/logging.conf')
logger = logging.getLogger('simulation_root')

readline.parse_and_bind('tab: complete')

parser = OptionParser()
parser.add_option("-f", "--file", dest="filename",
    help="write report to FILE", metavar="FILE")
parser.add_option("-p", "--plot",
    action="store_true", dest="create_plots", default=False,
    help="create some informative debugging plots")
options, args = parser.parse_args()
if options.filename is None:
    print 'Please submit a filename. "', sys.argv[0],'-h" for help'
    sys.exit(1)


#~ basename = os.path.split(options.filename)
basename = options.filename
print "basename",basename

##############################################################################
############################# Actual program start ###########################
run = readcorsika.readCorsikaFile( options.filename )

################ run = absorb_photons.absorb_photons(run)#####################    
spw = fact_si.interpolate_from_files( 
                fact_si.Wavelength_dependency_description_files )
data = run['photonData']
data = absorb_photons._absorb_photons( data, absorption_key='wavelength', absorption_func=spw)
data = turnIntoZenit.turn_photons_into_zenith( data )

isInInnerCircle = absorb_photons.IsInInnerCircle(210)
data = absorb_photons._absorb_photons( data, absorption_key='ground_pos_turned', absorption_func=isInInnerCircle)


#run = {}
#run['photonData'] = {}
#pos = np.zeros((10000,3))
#dir = np.zeros((10000,3))

#dir[:,2] = 1.
#pos[:,0] = np.random.uniform(-250,250,10000)
#pos[:,1] = np.random.uniform(-250,250,10000)

#run['photonData']['ground_pos_turned'] = pos
#run['photonData']['ground_dir_turned'] = dir



###################run = reflector.reflector(run)#############################
run = reflector.pushMirrorInfo( run )

run = reflector.findMirrorIntersection( run, 
        dictname='photonData', 
        pos_key='ground_pos_turned', 
        dir_key='ground_dir_turned',
        mirror_key='mirrors',
        new_pos_key='pos_on_mirror', 
        new_mirror_id_key='mirror_id', 
        new_time_key='time_to_mirror')

run = reflector.reflectOnMirrors( run, 
        dictname='photonData', 
        mirror_alignment_point_key='mirror_alignment_point',
        pos_key='pos_on_mirror', 
        dir_key='ground_dir_turned',
        new_angle_key='angle_on_mirror', 
        new_dir_key='dir_on_mirror' )
    
run = reflector.findFocalPlaneIntersection( run, 
        dictname='photonData', 
        focal_plane_key='focal_plane',
        pos_key='pos_on_mirror', 
        dir_key='dir_on_mirror',
        new_pos_key='pos_in_focal_plane', 
        new_time_key='time_to_focal_plane')
        
run = reflector.findPixelChid( run, 
        dictname='photonData', 
        focal_plane_key='focal_plane',
        pos_key='pos_in_focal_plane', 
        new_chid_key='chid', 
        new_Coordinator_key='Coordinator', 
        new_pixel_center_key='pixel_center')

if options.create_plots:
    readcorsika.make_plots(run)
    reflector.make_plots(run)

#import helpers as old_helpers
#old_helpers.plot_chids( run )

from detektorsim import FasterPulse
import detektorsim
import plotters
import matplotlib.pyplot as plt

pulse = FasterPulse()

pd = run['photonData']
event_headers=run['corsikaEventHeaders']
    
# find max event_id
event_ids = pd['event_id']
max_event_id = int(event_ids.max())


from fitswriter import RawDataFitsWriter, DrsCalibFitsWriter
drsfile = DrsCalibFitsWriter(basename+".drs.fits")
drsfile.write()
del drsfile

fitsfile = RawDataFitsWriter(basename+".fits")


#~ plt.ion()
#~ camplotter = plotters.CamPlotter('cam', dpi=80, s=80)
#~ fig2 = plt.figure()
PulseIntegral = detektorsim.measureIntegral()[0]

for event_id in range(max_event_id):
    #print "event_id",event_id
    ids = np.where(pd['event_id'] == event_id)[0]
    if not len(ids):
        continue
    chids = pd['chid'][ids].astype(int)
    times = pd['time'][ids]-pd['time_to_mirror'][ids]+pd['time_to_focal_plane'][ids]
    times -= times.min()
    ## The times are in ns, so at this point in convert them to sec.
    times *= 1e-9
    
    chidset = set(chids.tolist())
    
    timeLists = {}
    #pulseLists = {}
    for chid,t in zip(chids,times):
        if chid not in timeLists:
            timeLists[chid] = [t]
            #pulseLists[chid] = [Pulse(t)]
        else:
            timeLists[chid].append(t)
            #pulseLists[chid].append(Pulse(t))
    
    ########### Most simple Trigger ##########################################
    #~ maxNumPhotonsInPixel =0
    #~ for chid in timeLists:
        #~ if len(timeLists[chid]) > maxNumPhotonsInPixel:
            #~ maxNumPhotonsInPixel = len(timeLists[chid])
    #~ if not maxNumPhotonsInPixel > 4:
        #~ print "Discarding Event, because its boring"
        #~ continue
    #~ ######################## End of most simple trigger ######################
    
    # I try to use only SI units here, so one doesn't need to think...
    StartTimeOffset = -10 * 1e-9    # sec
    SimulatedTime = 512 * 1e-9      # sec
    SamplingFreqency = 2. * 1e9     # Hz
    
    NumberOfSamples = int(SimulatedTime*SamplingFreqency)
    StartTimeOffsetSamples = int(StartTimeOffset*SamplingFreqency)
    
    ########################## NSB and Dark noise addition ###################
    # For all Pixel
    NsbRate = 0.    # Mega Hz
    DcRate = 10.    # Mega Hz
    Rate = (NsbRate + DcRate) * 1e6 # in Hz now
    Period = 1./Rate                # in sec
    
    #print "adding NSB and Darknoise"
    for chid in range(1440):
        #print "chid:",chid, "event_id:",i
        # Generate Background Photons following an exponential with the 
        # expectation value *Period*

        # Expectation value of Photons in the simulated Time * 5 (contingency)
        numberOfNoisePhotons = int(SimulatedTime * Rate)
        timesOfNoisePhotons = np.random.exponential(Period, size=5*numberOfNoisePhotons)
        
        # Since there has to be some time simulated before the start of the Cherenkov Event,
        # I simply assume 1 Period of SimulatedTime before the arrival of the first Photon.
        timesOfNoisePhotons = timesOfNoisePhotons.cumsum() - SimulatedTime
        # timesOfNoisePhotons are now times in ns.
        
        # Now I cut away all Noise Photons which arrived after SimulatedTime
        timesOfNoisePhotons =timesOfNoisePhotons[timesOfNoisePhotons <= SimulatedTime]
        
        # The rest of the photons, I have in hand are actually falling into this pixel
        for t in timesOfNoisePhotons:
            if chid not in timeLists:
                timeLists[chid] = [t]
                #pulseLists[chid] = [Pulse(t)]
            else:
                timeLists[chid].append(t)
                #pulseLists[chid].append(Pulse(t))
    
    ##########################################################################
    
    
    
    rawdata = np.zeros( (1440, NumberOfSamples) )
    x = np.linspace(StartTimeOffset,SimulatedTime-StartTimeOffset,NumberOfSamples,False)
    #print "filling rawdata with pulses"
    for chid in timeLists:
        for t in timeLists[chid]:
            #print "adding noise: chid:",chid,"pulse:",pulse
            pulse(t,rawdata[chid])
    
    ##################################### AC coupling ########################
    
    ac_effect = PulseIntegral * numberOfNoisePhotons / SimulatedTime 
    #print "AC coupling effect :", ac_effect
    rawdata -= ac_effect
    
    
    ############### Electronics noise ########################################
    
    NoiseLevel = 0.25 # This is in termes of single photons
    rawdata += np.random.normal( loc=0., scale=NoiseLevel, size=rawdata.shape)
    
    
    ############ ADC simulation, i.e. digitization ###########################
    
    AdcFullScale = 200.      # Photons
    Adcoffset = 10.          # Photons
    AdcResolution = 4096.    # LSBs
    
    adcdata = rawdata.copy()
    adcdata += Adcoffset
    adcdata *= AdcResolution/AdcFullScale
    # move rawdata -2048 down, so it falls into the range of int16
    adcdata -= 2048
    adcdata = adcdata.astype('int16')
    
    
    
    
    
    ph_in_pix = np.zeros(1440)
    pixels = pd['chid'][ids].astype(int)
    for pix in pixels:
        ph_in_pix[pix] += 1
    #print "event=%d  photons=%d" % (event_id, len(ids))
    #print pixels
    #print ("energy=%f GeV" % event_headers[event_id]['total energy in GeV'])
    
    #calculating core distace:
    core = np.array(event_headers[event_id]['core location for scattered events in cm: (x,y)'][0])
    core_dist= np.sqrt((core*core).sum())
    #print "core dist=%f m" % (core_dist/100,)
    
    #~ camplotter(ph_in_pix)
    #~ plt.xlim(-22,22)
    #~ plt.ylim(-22,22)
    #~ 
    #~ plt.figure(fig2.number)
    #~ plt.cla()
    #~ for chid in chidset:
        #~ plt.plot(x,adcdata[chid],'.:',label=str(chid), hold=True)
    #~ plt.legend()
    #~ 
    
    
    
    ############### Save event to fitsfile ###################################
        # Format
        #
        #~ ('EventNum', '<u4'), 
        #~ ('TriggerNum', '<u4'), 
        #~ ('TriggerType', '<u4'), 
        #~ ('NumBoards', '<u4'), 
        #~ ('Errors', '|u1', 4), 
        #~ ('SoftTrig', '<u4'), 
        #~ ('UnixTimeUTC', '<u4', 2), 
        #~ ('BoardTime', '<u4', 40), 
        #~ ('StartCellData', '<u2', 1440), 
        #~ ('StartCellTimeMarker', '<u2', 160), 
        #~ ('Data', '<i2', (1440,1024))]) )
    
    fitsfile.rawdata['Data'] = adcdata
    fitsfile.rawdata['EventNum'] = event_id
    fitsfile.rawdata['NumBoards'] = 40
    
    fitsfile.write()
    

    #raw_input("?")
    
logging.info("the end.")    
    
    
    
    

